% function realRateCounterfactuals 
clear all; close all;
clear compareStru;
cucd=cd;


%% flagStru 
% .aZeroSS = 1 if assume a(0) is at SS 
%          = 0 otherwise 

%% 1. Load path and m-file 
paths.load='O:\PROJ_LIB\LJM\exog\BMJ\DtrFGSmoothInf3Split\drift noLabor 90 start 08q4 Pr7_2\end 2013.25 01 17 14'; 
%paths.load='O:\PROJ_LIB\LJM\exog\BMJ\DtrFGSmoothInf3Split\drift noLabor 90 start 08q4 Pr8\end 2013.25 01 21 14'; 

%paths.mfile='counterSettings';

histStru.extract={'pPCE'}; 
histStru.names={'CPI Core'};
histStru.scales=[4]; 

% Jumping point, with full state 
histStru.jumpDate=2008.5; 
% Date at which to decompose 
histStru.targetDate=2013.25; 

%% Plot settings 
figStru.startDate=1990.00; 
figStru.endDate  =2013.25; 

groupStruct.g1={'Inflation Drift'};
groupStruct.g2={'wCE me','wAHE me','pJGDP me','pCPI me','pPCE me','AvInf me'};
groupStruct.g3={'Labor Disutility','Price Markup','Wage Markup'};
groupStruct.g4={'G+NX+CIV','Risk'};
groupStruct.g5={'Neutral technology','MEI'};
groupStruct.g6={'Target Factor','Path Factor'}; 

histStru.groupnames={'Drift','Idiosyncratic','Markups',...
    'Demand','Tech','Policy','Residual'};

removeStru.A.groups={'Drift'}; 
removeStru.A.name  ='Detrended'; 
removeStru.A.nameComp='Inflation drift';   
removeStru.A.title ='Role of Inflation Drift'; 

removeStru.B.groups={'Drift','Idiosyncratic'}; 
removeStru.B.name='Detrended less ID'; 
removeStru.B.nameComp='Trend and ID'; 
removeStru.B.title ='Role of Idiosyncratic'; 

removeStru.C.groups={'Drift','Idiosyncratic','Residual'}; 
removeStru.C.name='Fundamental'; 
removeStru.C.nameComp='Trend, ID and Residual'; 
removeStru.C.title='Role of Residual'; 

%% Leftover.D = {'D+S+Policy'}; 
removeStru.D.groups={'Drift','Idiosyncratic','Residual','Markups'};  
removeStru.D.nameComp='Detrended less ID, Residual, Markup'; 
removeStru.D.name='Policy+Supply+Demand'; 
removeStru.D.title='Role of Markups'; 

%% Leftover.E={'S+Markups'}; 
removeStru.E.groups={'Drift','Idiosyncratic','Residual','Markups','Demand'}; 
removeStru.E.name='Only Policy and Supply'; 
removeStru.E.nameComp='Trend, ID, Residual and Demand'; 
removeStru.E.title='Role of Demand'; 

removeStru.F.groups={'Drift','Idiosyncratic','Residual','Markups','Demand','Policy'}; 
removeStru.F.name='Only Supply'; 
removeStru.F.title='Role of Policy'; 
removeStru.F.nameComp='Trend, ID, Residual, Demand and Policy'; 



%% 2. Load modelStructures 
cd(paths.load);
load modelStructures; 
%eval([paths.mfile]);
cd(cucd);

%% Unconditional Decomposition across Groups 
[uncondAll.counterMat,smoothStReport,dim]=...
    distillUnconditionalSub(KFStru,histStru.extract,names,groupStruct,histStru.groupnames); 

%% Removal of Groups 
dim.remove=length(fieldnames(removeStru)); 
leftOverMatrix=zeros(dim.obs,dim.extract,dim.remove+1); 

leftOverCell  =cell(dim.remove+1,1); 
leftOverCell(1)={'Smooth'};   


histStru.scaleMat=repmat( reshape( histStru.scales , [1 dim.extract] ) , dim.obs,  1 ); 
smoothStReport=histStru.scaleMat.*smoothStReport; 

removedMat=zeros(dim.obs,dim.extract,dim.remove); 
leftOverMat=zeros(dim.obs,dim.extract,dim.remove); 
removedCell=cell(dim.remove,1); 

removeStru.F.nameComp
tempFields=fieldnames(removeStru); 
for gInd=1:dim.remove; 
    fldnm=char(tempFields{gInd});
    removeStru.(fldnm).posRemoved=cellposition( removeStru.(fldnm).groups,...
        histStru.groupnames);
    removeStru.(fldnm).posLeftOver=setdiff( (1:dim.groups)' , removeStru.(fldnm).posRemoved );
    removeStru.(fldnm).removed=histStru.scaleMat.*( sum( ...
        uncondAll.counterMat(:,:,removeStru.(fldnm).posRemoved) ,3 ) );
    removeStru.(fldnm).leftover=smoothStReport-removeStru.(fldnm).removed;     
    removedMat(:,:,gInd) =removeStru.(fldnm).removed; 
    leftOverMat(:,:,gInd)=removeStru.(fldnm).leftover; 
    leftOverCell(gInd+1)={removeStru.(fldnm).name}; 
    removedCell(gInd)={removeStru.(fldnm).nameComp}; 
end 
clear temp*; 

leftOverMatrix(:,:,2:end)=leftOverMat; 
leftOverMatrix(:,:,1)=smoothStReport; 

cd(paths.load); 

for serind=1:dim.extract
    tempCell=crtcell( squeeze( leftOverMatrix(:,serind,:) ),...
        num2cprec( sampleStru.vector) , leftOverCell , [histStru.extract{serind},' leftover'] );
    xlswrite(['Counter ',strdate],tempCell,[histStru.extract{serind},' leftover'] );
    clear temp*
    
    tempCell=crtcell( squeeze( removedMat(:,serind,:) ),...
        num2cprec( sampleStru.vector) , removedCell , [histStru.extract{serind},' removed']);
    xlswrite(['Counter ',strdate],tempCell,[histStru.extract{serind},' removed'] );
    
end



%% Plotting Settings 
figStru.posStart=find(  sampleStru.vector == figStru.startDate); 
figStru.posEnd=find(  sampleStru.vector == figStru.endDate); 
figStru.posSample=(figStru.posStart:figStru.posEnd); 
figStru=plotAssignDefaults(figStru); 
figStru.xaxis =sampleStru.vector(figStru.posSample);
figStru.layout=[dim.extract 1]; 

% figStru.legend={'smooth','detrened (less drift)'}; 
% [~,~,~,figVec(1)]=...
%     distillUnconditionalRemove(KFStru,histStru.extract,names,groupStruct,histStru.groupnames,...
%     removeStru.A.groups,figStru); 


figVec=zeros(10,1); 


%% Plot Smooth vs. No Drift 
tempNameStru.states=histStru.names; 
tempNameStru.legend={'data',removeStru.A.name}; 
tempNameStru.title=cell(2,1); 
tempNameStru.title(1)={removeStru.A.title};
tempNameStru.title(2)={'Difference Solid vs. Dashed'};

figVec(1)=plot2StateMatsDiff(smoothStReport(figStru.posSample,:),...
    removeStru.A.leftover(figStru.posSample,:),...
    tempNameStru,figStru.xaxis ,figStru,'Detrending'); 


names.removeStru=fieldnames(removeStru); 
for rInd=2:dim.remove 
    fldnmCurr=names.removeStru{rInd}; 
    fldnmPast=names.removeStru{rInd-1}; 
    tempNameStru.legend={removeStru.(fldnmPast).name,...
        removeStru.(fldnmCurr).name};
    tempNameStru.title(1)={removeStru.(fldnmCurr).title};

    figVec(rInd)=plot2StateMatsDiff(removeStru.(fldnmPast).leftover(figStru.posSample,:),...
        removeStru.(fldnmCurr).leftover(figStru.posSample,:),...
        tempNameStru,figStru.xaxis ,figStru,'ID');

end 

tempNameStru.states=histStru.names; 
tempNameStru.legend={removeStru.C.name,removeStru.F.name}; 
tempNameStru.title={'Fundamental with Markups vs. Supply Only','Deviations from Supply'}; 
figVec(8)=plot2StateMatsDiff(removeStru.C.leftover(figStru.posSample,:),...
    removeStru.F.leftover(figStru.posSample,:),...
    tempNameStru,figStru.xaxis ,figStru,'Supply Only'); 

figVec=figVec( figVec~=0 ); 
cd(paths.load); 
delete sequential*; 
for ii=1:length(figVec); 
    print( figVec(ii),'-dpsc','sequential Decomp','-append'); 
end 
cd(cucd); 


dbstop; 

%% No Drift vs. No Drift ID
% removeStru.B.groups={'Drift','Idiosyncratic'}; 
% removeStru.B.name='Less Drift & ID'; 
tempNameStru.states=histStru.names; 


%% No Drift ID vs. No Drift Markups 
tempNameStru.states=histStru.names; 
tempNameStru.legend={removeStru.B.name,removeStru.C.name}; 
tempNameStru.title={'Role of "Residual"','Difference'}; 
figVec(3)=plot2StateMatsDiff(removeStru.B.leftover(figStru.posSample,:),...
    removeStru.C.leftover(figStru.posSample,:),...
    tempNameStru,figStru.xaxis ,figStru); 


%% No Drift-ID-Residual vs. No Drift-ID-Residual-
tempNameStru.states=histStru.names; 
tempNameStru.legend={removeStru.C.name,removeStru.D.name}; 
tempNameStru.title={'Left over is D+S+Policy','Leftover=markups'}; 
figVec(4)=plot2StateMatsDiff(removeStru.C.leftover(figStru.posSample,:),...
    removeStru.D.leftover(figStru.posSample,:),...
    tempNameStru,figStru.xaxis ,figStru,'Role of Markups'); 

tempNameStru.states=histStru.names; 
tempNameStru.legend={removeStru.D.name,removeStru.E.name}; 
tempNameStru.title={'Role of Demand','Difference=demand'}; 
figVec(5)=plot2StateMatsDiff(removeStru.D.leftover(figStru.posSample,:),...
    removeStru.E.leftover(figStru.posSample,:),...
    tempNameStru,figStru.xaxis ,figStru,'Supply and Policy Only'); 


tempNameStru.states=histStru.names; 
tempNameStru.legend={removeStru.E.name,removeStru.F.name}; 
tempNameStru.title={'Role of Policy','Difference=Supply Only'}; 
figVec(6)=plot2StateMatsDiff(removeStru.D.leftover(figStru.posSample,:),...
    removeStru.E.leftover(figStru.posSample,:),...
    tempNameStru,figStru.xaxis ,figStru,'Supply Only'); 



%% Plot States vs. Detrended 
figure; 
for ii=1:dim.extract; 
    subplot( graphRows, graphCols , ii);     
    tempMat=[smoothStReport(:,ii) removeStru.A.leftover(:,ii)]; 
    plot3samp(figStru.xaxis,tempMat(figStru.posSample,:),figStru); 
    title(char( histStru.names{ii} ) ); 
end 

for ii=1:dim.extract;    
        figVec(ii)=figure;    
        subplot(2,1,1);        
        figStru.flagNoLegend=0;
        %figStru.legend=names.legend;
        plot3samp(figStru.xaxis,[outStru.leftover(figStru.posSample,ii) ...
            stateMat(figStru.posSample,ii)],figStru);
        title('Left over vs. full');
        
        subplot(2,1,2);
        figStru.flagNoLegend=1;
        plot3samp(figStru.xaxis,...           
        outStru.removed(figStru.posSample,ii),...
            figStru);
        title('full-left over=removed');                 
        suptitle(char( extract{ii} ),14); 
end



plot2StateMats(removeStru.A.leftover(figStru.posSample,:),...
    removeStru.B.leftover(figStru.posSample,:),...
    tempNameStru,figStru.xaxis ,figStru,figStru.layout,'Detrended - idiosyncratic'); 










%% Remove the inflation drift 
unconStru.detrended=KFStru.smoothSt(:,histStru.pos)-squeeze( histStru.counterMat(:,:,1) ); 


%% Settings for conditional decomposition 
condStru.pos.eta=find( histStru.targetDate>= sampleStru.vector & sampleStru.vector > histStru.jumpDate); 
% pos.states: contains one more observation than eta, the initial jump
% state 
condStru.pos.states = find( histStru.targetDate>= sampleStru.vector & sampleStru.vector >=histStru.jumpDate);
condStru.jumpState=KFStru.smoothSt(find( sampleStru.vector == histStru.jumpDate    ) ,:);  
condStru.innovMat=KFStru.innovations(condStru.pos.eta,:); 
condStru.sample.forecast=sampleStru.vector( condStru.pos.eta );
condStru.sample.forWithInitial=sampleStru.vector( condStru.pos.states); 
condStru.Nforc =length(condStru.sample.forecast);
condStru.tauVec=model.addsol.tauVec(condStru.pos.eta); 


 distillForecastDecomp(model,condStru.tauVec,...
     condStru.Nforc,condStru.jumpState,condStru.innovMat,...
     histStru.extract,names,groupStruct,histStru.groupnames) 



%% 3. Compute baseline model 
%% Fields in momentBase now depends on number of splits 
%  .IRFStates   [ns nx  nper+1  nsplits]
%  .IRFObs      [nz nx  nper+1 nsplits] 
%  .corrObs.th  [nz nz ncor nsplits] 
%  .corrObs.sim [idem] (median)
%  .corrStates.th  [ns ns ncor nsplits] 
%  .corrStates.sim [idem] (median)
%  .STDObs.th   [nz nsplits] 
%  .STDObs.sim  [nz nsplits] (median)
%  .STDStates.th [ns nsplits]
%  .STDSstates.sim [ns splits] 
%  .varDecObs     [nz nx]
%  .varDecStates  [ns nx] 
[~,momentBase,KFBase,dim,spectrumBase]=modelAnalysisMultiple...
    (model.handle,model.param,model,sampleStru,compareStru); 

table.description=cell(4,2); 

%% 4. Assign aZero and Innovations
if flagStru.aZeroSS==true 
    counterStru.aZero=zeros( size(KFBase.aZero) ); 
else 
counterStru.aZero      =KFBase.aZero; 
end 
counterStru.innovations=KFBase.innovations; 
table.description(1,:)={'a(0) at SS?',flagStru.aZeroSS}; 

%% Additional dimensions 
names.variations=fieldnames(varStru); 
dim.variations=length(names.variations); 
dim.states=length(compareStru.stateDecomp); 
dim.IRFs=length(compareStru.stateIRF); 
dim.spectrum=length(compareStru.stateSpectrum);
dim.sample=length(sampleStru.vector); 
dim.moments=length(compareStru.stateMom); 

%% Add the function that will perform the counterfactuals 
model.addsol.counter=str2func([func2str( model.addsol.funcKfilter ),'Counterfactual']); 

compareStru.legend=cell(dim.variations+1,1); 
compareStru.legend(1)={'baseline'}; 

compareStru.states=zeros(dim.sample,dim.states,dim.variations+1);
compareStru.states(:,:,1)=KFBase.analysis.states; 
if dim.IRFs > 0
    %  .IRFStates   [ns nx  nper+1  nsplits]
    %  .IRFObs      [nz nx  nper+1 nsplits]
    compareStru.IRFObs=zeros( [dim.IRFObs dim.variations+1 ] );
    compareStru.IRFObs(:,:,:,1)=momentBase.IRFObs; 
    
    compareStru.IRFStates=zeros( [dim.IRFStates dim.variations+1 ] );
    compareStru.IRFStates(:,:,:,1)=momentBase.IRFStates;
end
if dim.moments > 0 
    compareStru.STDStates.th=zeros([dim.STDStates dim.variarions+1]); 
    compareStru.STDStates.sim=zeros([dim.STDStates dim.variarions+1]); 
  
    compareStru.STDStates.th(:,:,1)=momentBase.STDStates.th; 
    compareStru.STDStates.sim(:,:,1)=momentBase.STDStates.sim; 
end 
        
   
%%
for ii=1:dim.variations 
    fldnm=char(names.variations{ii});
    model.addsol=varStru.(fldnm).addsol;
    model.addsol.counter=str2func([func2str( model.addsol.funcKfilter ),'Counterfactual']);
    [~,varStru.(fldnm).moments,...
        varStru.(fldnm).KF,...
        varStru.(fldnm).spectrum,...
        varStru.(fldnm).counter]=...
        modelAnalysis(varStru.(fldnm).handle,varStru.(fldnm).param,...
        model,sampleStru,compareStru,counterStru);
    if isempty(varStru.(fldnm).KF)== false
        compareStru.states(:,:,ii+1)=varStru.(fldnm).KF.analysis.states;
    else
        compareStru.states(:,:,ii+1)=varStru.(fldnm).counter.states;        
    end
    compareStru.legend(ii+1)={char(varStru.(fldnm).name)}; 

    if dim.IRFs > 0
        compareStru.IRFs(:,:,:,ii+1)=varStru.(fldnm).moments.IRFStates;        
    end
    if dim.moments > 0 
        compareStru.moments.stdMed(:,ii+1)=...
            median( varStru.(fldnm).moments.STDSimStates, 2 ); 
    end 
end 

%% Extract desired states 
figStru=[];
figStru.start=3; 
figStru=plotAssignDefaults(figStru); 
figStru.legend=compareStru.legend;
figStru.xaxis =sampleStru.vector(figStru.start:end);
figStru.legendFontSize=14; 


paths.save=cr_dir(paths.load,compareStru.pathName);

states.figHandle=zeros(dim.states,1); 
for jj=1:dim.states; 
    states.figHandle(jj)=figure; 
    tempMat=squeeze( compareStru.states(figStru.start:end,jj,:) ); 
    plot3samp(figStru.xaxis,tempMat,figStru); 
    title(char( compareStru.stateDecomp{jj} ) ); 
    cd(paths.save); 
    
    table.states=crtcell(tempMat,num2cprec(figStru.xaxis),...
        compareStru.legend,compareStru.stateDecomp{jj}); 
    xlswrite('table',table.states,compareStru.stateDecomp{jj}); 

end 




if isempty(compareStru.pathName)==true 
    compareStru.pathName=['compare ',strdate]; 
else 
    compareStru.pathName=strcat(char(compareStru.pathName),' compare ',strdate); 
end 

cd(paths.save); 
delete statesFig*
for jj=1:dim.states 
print(states.figHandle(jj),'-dpsc','statesFig','-append');
end 
cd(cucd); 
disp('Done'); 


if compareStru.save.tables==true
    table.moments=crtcell(compareStru.moments.stdMed,compareStru.stateMom,...
        compareStru.legend,'Simulated Moments'); 

    printcell( table.moments ); 
    cd(paths.save); 
    xlswrite('table',table.moments,'moments'); 
end


IRFHandle=plot_irfperturb(compareStru.IRFs,[],[],...
    compareStru.legend,compareStru.stateIRF,names.shocks);
if compareStru.save.IRFs==true
    cd(paths.save);
    delete IRF*;
    for jj=1:size(compareStru.IRFs,2)
        print(IRFHandle(jj),'-dpsc','IRFs','-append');
    end
    cd(cucd);
    
end